Root_Atlas.rds (or get it by running through notebook 2, 3, 4, 5, 6, 7 & 8)
Ground_Tissue_Atlas.rds (get it by running through notebook 7)
Ground_Tissue_Atlas_latent_time.txt (get it by running through notebook 8)
Epidermis_LRC_Atlas.rds (get it by running through notebook 7)
Epidermis_LRC_Atlas_latent_time.txt (get it by running through notebook 8, need to change tissue/lineage)
Columella_Atlas.rds (get it by running through notebook 7)
Columella_Atlas_latent_time.rds (get it by running through notebook 8, need to change tissue/lineage)
Stele_Atlas.rds (get it by running through notebook 7)
Stele_Atlas_latent_time.rds (get it by running through notebook 8, need to change tissue/lineage)
rm(list=ls())
# Set the working directory to where folders named after the samples are located.
# The folder contains spliced.mtx, unspliced.mtx, barcodes and gene id files, and json files produced by scKB that documents the sequencing stats.
setwd("/scratch/AG_Ohler/CheWei/scKB")
# Load libraries
suppressMessages(library(Seurat))
suppressMessages(library(CytoTRACE))
suppressMessages(library(RColorBrewer))
suppressMessages(library(ggplot2))
suppressMessages(library(grid))
sessionInfo()
# Plotting function for time
plot_time <- function(rc.integrated){
p1 <- FeaturePlot(end.cor.integrated, features = "CytoTRACE", pt.size=0.5)+ scale_colour_gradientn(colours = rev(brewer.pal(11,"Spectral")))
p2 <- FeaturePlot(end.cor.integrated, features = "latenttime", pt.size=0.5)+ scale_colour_gradientn(colours = brewer.pal(11,"Spectral"))+ggtitle("Latent Time")
p3 <- FeaturePlot(end.cor.integrated, features = "consensus.time", pt.size=0.5)+ scale_colour_gradientn(colours = brewer.pal(11,"Spectral"))+ggtitle("Consensus Time")
options(repr.plot.width=16, repr.plot.height=16)
gl <- lapply(list(p1, p2, p3), ggplotGrob)
gwidth <- do.call(unit.pmax, lapply(gl, "[[", "widths"))
gl <- lapply(gl, "[[<-", "widths", value = gwidth)
gridExtra::grid.arrange(grobs=gl, ncol=2)
}
# Read in atlas
rc.integrated <- readRDS('./Root_Atlas.rds')
# Read in ground tissue atlas
end.cor.integrated <- readRDS('./supp_data/Ground_Tissue_Atlas.rds')
# Function for normalizing values to range 0 to 1
range01 <- function(x) {
(x - min(x))/(max(x) - min(x))
}
# Read in latent time
latent.time <- as.numeric(read.table("./supp_data/Ground_Tissue_Atlas_latent_time.txt", header=F)$V1)
end.cor.integrated$latenttime <- latent.time
# Calculate and store consensus time
consensus.time <- range01((-(end.cor.integrated$CytoTRACE-1)+end.cor.integrated$latenttime)/2)
end.cor.integrated$consensus.time <- consensus.time
plot_time(end.cor.integrated)
# Correlation among CytoTRACE, latent time and consensus time
res <- cor(data.frame(CytoTRACE=end.cor.integrated$CytoTRACE,latent.time=end.cor.integrated$latenttime,consensus.time=end.cor.integrated$consensus.time))
round(res, 2)
# Save Seurat object
saveRDS(end.cor.integrated,'./supp_data/Ground_Tissue_Atlas.rds')
end.cor.integrated <- readRDS('./supp_data/Epidermis_LRC_Atlas.rds')
latent.time <- as.numeric(read.table("./supp_data/Epidermis_LRC_Atlas_latent_time.txt", header=F)$V1)
end.cor.integrated$latenttime <- latent.time
consensus.time <- range01((-(end.cor.integrated$CytoTRACE-1)+end.cor.integrated$latenttime)/2)
end.cor.integrated$consensus.time <- consensus.time
plot_time(end.cor.integrated)
res <- cor(data.frame(CytoTRACE=end.cor.integrated$CytoTRACE,latent.time=end.cor.integrated$latenttime,consensus.time=end.cor.integrated$consensus.time))
round(res, 2)
saveRDS(end.cor.integrated,'./supp_data/Epidermis_LRC_Atlas.rds')
end.cor.integrated <- readRDS('./supp_data/Columella_Atlas.rds')
latent.time <- as.numeric(read.table("./supp_data/Columella_Atlas_latent_time", header=F)$V1)
end.cor.integrated$latenttime <- latent.time
consensus.time <- range01((-(end.cor.integrated$CytoTRACE-1)+end.cor.integrated$latenttime)/2)
end.cor.integrated$consensus.time <- consensus.time
plot_time(end.cor.integrated)
res <- cor(data.frame(CytoTRACE=end.cor.integrated$CytoTRACE,latent.time=end.cor.integrated$latenttime,consensus.time=end.cor.integrated$consensus.time))
round(res, 2)
saveRDS(end.cor.integrated,'./supp_data/Columella_Atlas.rds')
end.cor.integrated <- readRDS('./supp_data/Stele_Atlas.rds')
latent.time <- as.numeric(read.table("./supp_data/Stele_Atlas_latent_time.txt", header=F)$V1)
end.cor.integrated$latenttime <- latent.time
end.cor.integrated$CytoTRACE[which(is.na(end.cor.integrated$CytoTRACE))] <- end.cor.integrated$latenttime[which(is.na(end.cor.integrated$CytoTRACE))]
consensus.time <- range01((-(end.cor.integrated$CytoTRACE-1)+end.cor.integrated$latenttime)/2)
end.cor.integrated$consensus.time <- consensus.time
plot_time(end.cor.integrated)
res <- cor(data.frame(CytoTRACE=end.cor.integrated$CytoTRACE,latent.time=end.cor.integrated$latenttime,consensus.time=end.cor.integrated$consensus.time))
round(res, 2)
saveRDS(end.cor.integrated,'./supp_data/Stele_Atlas.rds')
# Load 4 tissue/lineages
end.cor.integrated <- readRDS('./supp_data/Ground_Tissue_Atlas.rds')
epi.integrated <- readRDS('./supp_data/Epidermis_LRC_Atlas.rds')
col.integrated <- readRDS('./supp_data/Columella_Atlas.rds')
stl.integrated <- readRDS('./supp_data/Stele_Atlas.rds')
# Cell index for each tissue/lineages
end.cor.traj.idx <- which(rc.integrated$celltype.anno == "Endodermis" | rc.integrated$celltype.anno == "Cortex" | rc.integrated$celltype.anno == "Putative Quiescent Center"| rc.integrated$celltype.anno == "Stem Cell Niche")
epi.traj.idx <- which(rc.integrated$celltype.anno == "Atrichoblast" | rc.integrated$celltype.anno == "Trichoblast" | rc.integrated$celltype.anno == "Lateral Root Cap" | rc.integrated$celltype.anno == "Putative Quiescent Center"| rc.integrated$celltype.anno == "Stem Cell Niche")
col.traj.idx <- which(rc.integrated$celltype.anno == "Columella" | rc.integrated$celltype.anno == "Putative Quiescent Center"| rc.integrated$celltype.anno == "Stem Cell Niche")
stl.traj.idx <- which(rc.integrated$celltype.anno == "Pericycle" | rc.integrated$celltype.anno == "Phloem" | rc.integrated$celltype.anno == "Xylem" | rc.integrated$celltype.anno == "Procambium" | rc.integrated$celltype.anno == "Putative Quiescent Center"| rc.integrated$celltype.anno == "Stem Cell Niche")
scn.traj.idx <- which(rc.integrated$celltype.anno == "Putative Quiescent Center"| rc.integrated$celltype.anno == "Stem Cell Niche")
# Combine consensus time
t1 <- rep(0,ncol(rc.integrated))
t2 <- rep(0,ncol(rc.integrated))
t3 <- rep(0,ncol(rc.integrated))
t4 <- rep(0,ncol(rc.integrated))
t1[end.cor.traj.idx] <- end.cor.integrated$consensus.time
t2[epi.traj.idx] <- epi.integrated$consensus.time
t3[col.traj.idx] <- col.integrated$consensus.time
t4[stl.traj.idx] <- stl.integrated$consensus.time
consensus.time <- t1+t2+t3+t4
# Take average of SCN and QC cells
consensus.time[scn.traj.idx] <- consensus.time[scn.traj.idx]/4
# Store final results of consensus time
rc.integrated$consensus.time <- consensus.time
# Plot consensus time
options(repr.plot.width=10, repr.plot.height=10)
FeaturePlot(rc.integrated, features = "consensus.time", pt.size=0.5)+ scale_colour_gradientn(colours = brewer.pal(11,"Spectral"))+ggtitle("Consensus Time")
# Divide cells into ten groups of even size based on consensus time
d <- round(ncol(rc.integrated)/10);
d
names(consensus.time) <- seq(1,ncol(rc.integrated),1)
consensus.time.group <- rep("Unknown",ncol(rc.integrated))
consensus.time.group[as.numeric(names(sort(consensus.time)[1:d]))] = "T0"
consensus.time.group[as.numeric(names(sort(consensus.time)[(d+1):(2*d)]))] = "T1"
consensus.time.group[as.numeric(names(sort(consensus.time)[(2*d+1):(3*d)]))] = "T2"
consensus.time.group[as.numeric(names(sort(consensus.time)[(3*d+1):(4*d)]))] = "T3"
consensus.time.group[as.numeric(names(sort(consensus.time)[(4*d+1):(5*d)]))] = "T4"
consensus.time.group[as.numeric(names(sort(consensus.time)[(5*d+1):(6*d)]))] = "T5"
consensus.time.group[as.numeric(names(sort(consensus.time)[(6*d+1):(7*d)]))] = "T6"
consensus.time.group[as.numeric(names(sort(consensus.time)[(7*d+1):(8*d)]))] = "T7"
consensus.time.group[as.numeric(names(sort(consensus.time)[(8*d+1):(9*d)]))] = "T8"
consensus.time.group[as.numeric(names(sort(consensus.time)[(9*d+1):ncol(rc.integrated)]))] = "T9"
table(consensus.time.group)
# Store consensus time group
rc.integrated$consensus.time.group <- consensus.time.group
# Plot consensus time group
options(repr.plot.width=10, repr.plot.height=10)
DimPlot(rc.integrated, reduction = "umap", group.by = "consensus.time.group", pt.size=0.5, cols=brewer.pal(10,"Spectral"))
saveRDS(rc.integrated,'./Root_Atlas.rds')
end.cor.integrated$consensus.time.group <- rc.integrated$consensus.time.group[end.cor.traj.idx];
epi.integrated$consensus.time.group <- rc.integrated$consensus.time.group[epi.traj.idx];
col.integrated$consensus.time.group <- rc.integrated$consensus.time.group[col.traj.idx];
stl.integrated$consensus.time.group <- rc.integrated$consensus.time.group[stl.traj.idx];
# Plot consensus time group for ground tissue
options(repr.plot.width=10, repr.plot.height=10)
DimPlot(end.cor.integrated, reduction = "umap", group.by = "consensus.time.group", pt.size=0.5, cols=brewer.pal(10,"Spectral"))
# Plot consensus time group for epidermis and LRC
options(repr.plot.width=10, repr.plot.height=10)
DimPlot(epi.integrated, reduction = "umap", group.by = "consensus.time.group", pt.size=0.5, cols=brewer.pal(10,"Spectral"))
# Plot consensus time group for columella
options(repr.plot.width=10, repr.plot.height=10)
DimPlot(col.integrated, reduction = "umap", group.by = "consensus.time.group", pt.size=0.5, cols=brewer.pal(10,"Spectral"))
# Plot consensus time group for stele
options(repr.plot.width=10, repr.plot.height=10)
DimPlot(stl.integrated, reduction = "umap", group.by = "consensus.time.group", pt.size=0.5, cols=brewer.pal(10,"Spectral"))
saveRDS(end.cor.integrated, './supp_data/Ground_Tissue_Atlas.rds')
saveRDS(epi.integrated, './supp_data/Epidermis_LRC_Atlas.rds')
saveRDS(col.integrated, './supp_data/Columella_Atlas.rds')
saveRDS(stl.integrated, './supp_data/Stele_Atlas.rds')